﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace LPLib
{
    /// <summary>
    /// Data structure to contain gapscore components.
    /// </summary>
    public class GapScores
    {
        /// <summary>
        /// The score for opening a new deletion.
        /// </summary>
        public double OpenDeletionScore = -3.912; // log(0.02), 0.02 = GapProbabilities.OpenDeletionProbability default value.
        /// <summary>
        /// The score for continuing a deletion for one more base.
        /// </summary>
        public double ContinueDeletionScore = -0.1054; // log(0.9), 0.9 = GapProbabilities.ContinueDeletionProbability default value.
        /// <summary>
        /// The score for opening a new insertion.
        /// </summary>
        public double OpenInsertionScore = -3.912;
        /// <summary>
        /// The score for continuing an insertion for one more base.
        /// </summary>
        public double ContinueInsertionScore = -0.1054;

        public GapScores() { } // uses default values

        public GapScores(GapProbabilities gapProbabilities)
        {
            OpenDeletionScore = Math.Log(gapProbabilities.OpenDeletionProbability);
            OpenInsertionScore = Math.Log(gapProbabilities.OpenInsertionProbability);
            ContinueDeletionScore = Math.Log(gapProbabilities.ContinueDeletionProbability);
            ContinueInsertionScore = Math.Log(gapProbabilities.ContinueInsertionProbability);
        }
    }

    public class GapProbabilities
    {
        /// <summary>
        /// The probability for opening a new deletion.
        /// </summary>
        public double OpenDeletionProbability = 0.02;
        /// <summary>
        /// The probability for continuing a deletion for one more base.
        /// </summary>
        public double ContinueDeletionProbability = 0.9;
        /// <summary>
        /// The probability for opening a new insertion.
        /// </summary>
        public double OpenInsertionProbability = 0.02;
        /// <summary>
        /// The probability for continuing an insertion for one more base.
        /// </summary>
        public double ContinueInsertionProbability = 0.9;
    }

    /// <summary>
    /// Abstract class that specifies what methods alignment kernel must provide. It may compare Polymers of different nucleotide encoding types.
    /// </summary>
    /// <typeparam name="T">The type of the first nucleotide encoding.</typeparam>
    /// <typeparam name="U">The type of the second nucleotide encoding.</typeparam>
    public abstract class AlignmentKernelAsymmetric<T, U>
    {
        /// <summary>
        /// Reports the status of the symbol as a gap.
        /// </summary>
        /// <param name="n">The symbol to be tested.</param>
        /// <returns>True if the symbol represents a gap.</returns>
        abstract public bool IsGap(U n);
        /// <summary>
        /// Computes the comparison score for the alignment of the nucleotide at a given position in the query with the 
        /// nucleotide at a given position in the reference.
        /// </summary>
        /// <param name="posInQuery">The position in query.</param>
        /// <param name="posInReference">The position in the reference.</param>
        /// <returns></returns>
        abstract public double ComparisonScore(int posInQuery, int posInReference);
        /// <summary>
        /// Computes the score for opening an insertion in the query at a given query position when that position 
        /// is aligned against a given position in the reference.
        /// </summary>
        /// <param name="posInQuery">The query position.</param>
        /// <param name="posInReference">The reference position.</param>
        /// <returns>The cost of opening an insertion at the given coordinates.</returns>
        abstract public double OpenInsertionScore(int posInQuery, int posInReference);
        /// <summary>
        /// Computes the score for continuing an insertion in the query at a given query position when that position 
        /// is aligned against a given position in the reference.
        /// </summary>
        /// <param name="posInQuery">The query position.</param>
        /// <param name="posInReference">The reference position.</param>
        /// <returns>The cost of opening an insertion at the given coordinates.</returns>
        abstract public double ContinueInsertionScore(int posInQuery, int posInReference);
        /// <summary>
        /// Computes the score for opening a deletion in the query at a given query position when that position 
        /// is aligned against a given position in the reference.
        /// </summary>
        /// <param name="posInQuery">The query position.</param>
        /// <param name="posInReference">The reference position.</param>
        /// <returns>The cost of opening an insertion at the given coordinates.</returns>
        abstract public double OpenDeletionScore(int posInQuery, int posInReference);
        /// <summary>
        /// Computes the score for continuing a deletion in the query at a given query position when that position 
        /// is aligned against a given position in the reference.
        /// </summary>
        /// <param name="posInQuery">The query position.</param>
        /// <param name="posInReference">The reference position.</param>
        /// <returns>The cost of opening an insertion at the given coordinates.</returns>
        abstract public double ContinueDeletionScore(int posInQuery, int posInReference);
        /// <summary>
        /// Sets the query Polymer.
        /// </summary>
        /// <param name="query">The Polymer to set as query.</param>
        abstract public void SetQuery(Polymer<T> query);
        /// <summary>
        /// Gets the query.
        /// </summary>
        /// <returns>The query Polymer.</returns>
        abstract public Polymer<T> GetQuery();
        /// <summary>
        /// Sets the reference Polymer.
        /// </summary>
        /// <returns>The Polymer to set as reference.</returns>
        abstract public void SetReference(Polymer<U> reference);
        /// <summary>
        /// Gets the reference.
        /// </summary>
        /// <returns>The reference Polymer.</returns>
        abstract public Polymer<U> GetReference();
        /// <summary>
        /// The score taken as the default for the query at the given position.
        /// </summary>
        /// <param name="posInQuery">The position in the query.</param>
        /// <returns>The default score.</returns>
        abstract public double PriorQueryScore(int posInQuery);
        /// <summary>
        /// The score taken as default for the query.
        /// </summary>
        /// <returns>The default query score.</returns>
        abstract public double PriorQueryScore();
        /// <summary>
        /// Gets the array of default query scores.
        /// </summary>
        /// <returns>An array of doubles giving the default query score by position.</returns>
        abstract public double[] PriorQueryScoreArray();
    }

    public class MNMNAlignmentKernelAsymmetric : AlignmentKernelAsymmetric<byte, byte>
    {
        private NucleotideEvolutionModel model;

        private Polymer<byte> query;
        private Polymer<byte> reference;
        private GapScores gapScores;
        private double distance;
        public double Distance
        {
            get { return distance; }
            set
            {
                distance = value;
                model.Time = distance;
                computeComparisonScoreArray();
                mutationCost = Math.Log(1 - Math.Exp(-distance));
            }
        }

        private double mutationCost;
        public GapScores GapScores
        {
            set { gapScores = value; }
        }

        private double[,] comparisonScore;

        public MNMNAlignmentKernelAsymmetric(NucleotideEvolutionModel model)
        {
            this.model = model;
            gapScores = new GapScores();
            this.distance = model.Time;
            mutationCost = Math.Log(1 - Math.Exp(-distance));
            computeComparisonScoreArray();
        }

        public MNMNAlignmentKernelAsymmetric(NucleotideEvolutionModel model, Polymer<byte> query, Polymer<byte> reference)
        {
            this.query = query;
            this.reference = reference;
            this.GapScores = new GapScores();
            this.distance = model.Time;
            mutationCost = Math.Log(1 - Math.Exp(-distance));
            computeComparisonScoreArray();
        }

        public MNMNAlignmentKernelAsymmetric(NucleotideEvolutionModel model, PolymerPair<byte, byte> pair)
        {
            this.query = pair.Polymer0;
            this.reference = pair.Polymer1;
            this.gapScores = new GapScores();
            this.distance = model.Time;
        }

        private void computeComparisonScoreArray()
        {
            comparisonScore = new double[MixedNucleotideByte.NSymbols, MixedNucleotideByte.NSymbols];
            for (byte bR = 0; bR < MixedNucleotideByte.NSymbols; bR++)
            {
                double[] pR = TypeConverter.ToPMF(bR);
                for (byte bQ = 0; bQ < MixedNucleotideByte.NSymbols; bQ++)
                {
                    double[] pQ = TypeConverter.ToPMF(bQ);

                    for (int r = 0; r < 4; r++) // pure nucleotide states in reference sequence
                    {
                        double x = 0;
                        for (int q = 0; q < 4; q++) // pure nuclelotide states in query sequence
                        {
                            x += pQ[q] * model.TransitionProbability(r, q);
                        }
                        if (x > 0)
                        {
                            comparisonScore[bQ, bR] += pR[r] * Math.Log(x);
                        }
                        else if (pR[r] > 0)
                        {
                            comparisonScore[bQ, bR] += double.NegativeInfinity;
                        }
                        // else add zero = do nothing.
                    }
                }
            }
        }

        override public bool IsGap(byte n)
        {
            return n == 16;
        }

        override public double ComparisonScore(int posInQuery, int posInReference)
        {
            return comparisonScore[query.Seq[posInQuery], reference.Seq[posInReference]];
        }

        override public double OpenInsertionScore(int posInQuery, int posInReference)
        {
            return mutationCost + gapScores.OpenInsertionScore + PureNucleotideInt.LogPrior;
        }

        override public double ContinueInsertionScore(int posInQuery, int posInReference)
        {
            return gapScores.ContinueInsertionScore + PureNucleotideInt.LogPrior;
        }

        override public double OpenDeletionScore(int posInQuery, int posInReference)
        {
            return mutationCost + gapScores.OpenDeletionScore;
        }

        override public double ContinueDeletionScore(int posInQuery, int posInReference)
        {
            return gapScores.ContinueDeletionScore;
        }

        override public void SetQuery(Polymer<byte> query)
        {
            this.query = query;
        }

        public override Polymer<byte> GetQuery()
        {
            return query;
        }
        override public void SetReference(Polymer<byte> reference)
        {
            this.reference = reference;
        }

        public override Polymer<byte> GetReference()
        {
            return reference;
        }

        override public double PriorQueryScore(int posInQuery)
        {
            return PureNucleotideInt.LogPrior;
        }

        override public double PriorQueryScore()
        {
            return PureNucleotideInt.LogPrior;
        }

        override public double[] PriorQueryScoreArray()
        {
            double[] returnValue = new double[query.Seq.Length];

            returnValue[0] = PureNucleotideInt.LogPrior;
            for (int i = 1; i < query.Seq.Length; i++)
            {
                returnValue[i] = PureNucleotideInt.LogPrior + returnValue[i - 1];
            }
            return returnValue;
        }
    }

    public class PMFAlignmentKernel : AlignmentKernelAsymmetric<double[], double[]>
    {
        private Polymer<double[]> query;
        private Polymer<double[]> reference;
        private GapScores gapScores;
        private double mutationCost = MixedNucleotideByte.LogPrior;

        public GapScores GapScores
        {
            set { gapScores = value; }
        }

        public PMFAlignmentKernel(GapScores gapScores)
        {
            this.gapScores = gapScores;
        }

        public PMFAlignmentKernel(Polymer<double[]> query, Polymer<double[]> reference)
        {
            this.query = query;
            this.reference = reference;
            this.gapScores = new GapScores();
        }

        public PMFAlignmentKernel(PolymerPair<double[], double[]> pair)
        {
            this.query = pair.Polymer0;
            this.reference = pair.Polymer1;
            this.gapScores = new GapScores();
        }

        override public bool IsGap(double[] n)
        {
            return n[4] == 1;
        }

        override public double ComparisonScore(int posInQuery, int posInReference)
        {
            return MixedNucleotideByte.LogPrior + Math.Log(NucleotidePMF.IDFlux(query.Seq[posInQuery], reference.Seq[posInReference]));
        }

        override public double OpenInsertionScore(int posInQuery, int posInReference)
        {
            return mutationCost + gapScores.OpenInsertionScore + 2 * PureNucleotideInt.LogPrior;
        }

        override public double ContinueInsertionScore(int posInQuery, int posInReference)
        {
            return gapScores.ContinueInsertionScore + 2 * PureNucleotideInt.LogPrior;
        }

        override public double OpenDeletionScore(int posInQuery, int posInReference)
        {
            return mutationCost + gapScores.OpenDeletionScore + PureNucleotideInt.LogPrior;
        }

        override public double ContinueDeletionScore(int posInQuery, int posInReference)
        {
            return gapScores.ContinueDeletionScore + PureNucleotideInt.LogPrior;
        }

        override public void SetQuery(Polymer<double[]> query)
        {
            this.query = query;
        }

        public override Polymer<double[]> GetQuery()
        {
            return query;
        }
        override public void SetReference(Polymer<double[]> reference)
        {
            this.reference = reference;
        }

        public override Polymer<double[]> GetReference()
        {
            return reference;
        }

        override public double PriorQueryScore(int posInQuery)
        {
            return PureNucleotideInt.LogPrior;
        }

        override public double PriorQueryScore()
        {
            return PureNucleotideInt.LogPrior;
        }

        override public double[] PriorQueryScoreArray()
        {
            double[] returnValue = new double[query.Seq.Length];

            returnValue[0] = PureNucleotideInt.LogPrior;
            for (int i = 1; i < query.Seq.Length; i++)
            {
                returnValue[i] = PureNucleotideInt.LogPrior + returnValue[i - 1];
            }
            return returnValue;
        }
    }

    public class LikelihoodByteAlignmentKernelAsymmetric : AlignmentKernelAsymmetric<double[], byte>
    {
        private Polymer<double[]> query;
        private Polymer<byte> reference;
        private GapScores gapScores;
        private PolymerPair<double[], byte> pair;
        private double[][] loglikelihood;
        private double distance;

        public double Distance
        {
            get { return distance; }
            set
            {
                distance = value;
                mutationCost = Math.Log(1 - Math.Exp(-distance));
            }
        }

        private double mutationCost;

        public GapScores GapScores
        {
            set { gapScores = value; }
        }

        public LikelihoodByteAlignmentKernelAsymmetric(double distance)
        {
            this.distance = distance;
            mutationCost = Math.Log(1 - Math.Exp(-distance));
            this.gapScores = new GapScores();
            pair = new PolymerPair<double[], byte>();
        }

        public LikelihoodByteAlignmentKernelAsymmetric(Polymer<double[]> query, Polymer<byte> reference, double distance)
        {
            this.distance = distance;
            mutationCost = Math.Log(1 - Math.Exp(-distance));
            this.query = query;
            this.reference = reference;
            this.gapScores = new GapScores();
            makeLogLikelihood();
        }

        public LikelihoodByteAlignmentKernelAsymmetric(PolymerPair<double[], byte> pair, double distance)
        {
            this.distance = distance;
            mutationCost = Math.Log(1 - Math.Exp(-distance));
            this.query = pair.Polymer0;
            this.reference = pair.Polymer1;
            this.gapScores = new GapScores();
            makeLogLikelihood();
        }

        private void makeLogLikelihood()
        {
            double[] pb = new double[5];
            loglikelihood = new double[query.Seq.Length][];
            for (int pos = 0; pos < query.Seq.Length; pos++)
            {
                loglikelihood[pos] = new double[17];
                for (byte b = 0; b < 16; b++) // reference degenerate nucleotides
                {
                    pb = TypeConverter.ToPMF(b);
                    double averageloglikelihood = 0;
                    for (int iNuc = 0; iNuc < 5; iNuc++)
                    {
                        if (query.Seq[pos][iNuc] > 0)
                        {
                            averageloglikelihood += pb[iNuc] * Math.Log(query.Seq[pos][iNuc]);
                        }
                    }
                    loglikelihood[pos][b] = averageloglikelihood;
                }
                loglikelihood[pos][16] = Math.Log(query.Seq[pos][4]);
            }
        }

        override public bool IsGap(byte n)
        {
            return n == 16;
        }

        override public double ComparisonScore(int posInQuery, int posInReference)
        {
            return loglikelihood[posInQuery][reference.Seq[posInReference]];
        }

        override public double OpenInsertionScore(int posInQuery, int posInReference)
        {
            return loglikelihood[posInQuery][16];
            // compare: mutationCost + openInsertionScore + logPrior
        }

        override public double ContinueInsertionScore(int posInQuery, int posInReference)
        {
            return loglikelihood[posInQuery][16] + gapScores.ContinueInsertionScore - gapScores.OpenInsertionScore;
            // compare: continueInsertionScore + logPrior
        }

        override public double OpenDeletionScore(int posInQuery, int posInReference)
        {
            return mutationCost + gapScores.OpenDeletionScore;
            // compare: same. But make sure mutationCost is equivalent
        }

        override public double ContinueDeletionScore(int posInQuery, int posInReference)
        {
            return gapScores.ContinueDeletionScore;
            // compare: same
        }

        override public void SetQuery(Polymer<double[]> query)
        {
            this.query = query;
            makeLogLikelihood();
        }

        public override Polymer<double[]> GetQuery()
        {
            return query;
        }

        override public void SetReference(Polymer<byte> reference)
        {
            this.reference = reference;
        }

        public override Polymer<byte> GetReference()
        {
            return reference;
        }

        override public double PriorQueryScore(int posInQuery)
        {
            return loglikelihood[posInQuery][15] + MixedNucleotideByte.LogPrior;
        }

        override public double PriorQueryScore()
        {
            return double.NaN;
        }

        override public double[] PriorQueryScoreArray()
        {
            double[] returnValue = new double[query.Seq.Length];

            returnValue[0] = PriorQueryScore(0);
            for (int i = 1; i < query.Seq.Length; i++)
            {
                returnValue[i] = PriorQueryScore(i) + returnValue[i - 1];
            }
            return returnValue;
        }

        public static double[] NScoreArray(Polymer<double[]> queryLikelihood)
        {
            double[] loglikelihood = new double[queryLikelihood.Seq.Length];
            double[] pmfForN = new double[5] { 0.25, 0.25, 0.25, 0.25, 0 };
            for (int pos = 0; pos < queryLikelihood.Seq.Length; pos++)
            {
                double likelihood = 0;
                likelihood = 0;
                for (int iNuc = 0; iNuc < 4; iNuc++)
                {
                    likelihood += queryLikelihood.Seq[pos][iNuc];
                }
                loglikelihood[pos] = Math.Log(likelihood) + PureNucleotideInt.LogPrior;
            }

            // make cumulative
            for (int pos = 1; pos < queryLikelihood.Seq.Length; pos++)
            {
                loglikelihood[pos] += loglikelihood[pos - 1];
            }
            return loglikelihood;
        }
    }

    //public class LikelihoodPMFAlignmentKernel : AlignmentKernel<double[], double[]>
    //{
    //    // for this comparer, the query is a likelihood array: the values at each position need not sum to one.
    //    private Polymer<double[]> query;
    //    private Polymer<double[]> reference;
    //    private GapScores gapScores;
    //    private PolymerPair<double[], double[]> pair;
    //    private double[][] loglikelihood;
    //    private double[] logSumLikelihood;

    //    private double distance;

    //    public double Distance
    //    {
    //        get { return distance; }
    //        set
    //        {
    //            distance = value;
    //            mutationCost = Math.Log(distance);
    //        }
    //    }
    //    private double mutationCost;

    //    public GapScores GapScores
    //    {
    //        set { gapScores = value; }
    //    }

    //    public LikelihoodPMFAlignmentKernel(double distance, GapScores gapScores)
    //    {
    //        this.distance = distance;
    //        this.mutationCost = Math.Log(distance);
    //        this.gapScores = gapScores;
    //        pair = new PolymerPair<double[], double[]>();
    //    }

    //    public LikelihoodPMFAlignmentKernel(Polymer<double[]> query, Polymer<double[]> reference, double distance, GapScores gapScores)
    //    {
    //        this.distance = distance;
    //        this.mutationCost = Math.Log(distance);
    //        this.query = query;
    //        this.reference = reference;
    //        this.gapScores = gapScores;
    //        makeLogLikelihood();
    //    }

    //    public LikelihoodPMFAlignmentKernel(PolymerPair<double[], double[]> pair, double distance, GapScores gapScores)
    //    {
    //        this.distance = distance;
    //        this.mutationCost = Math.Log(distance);
    //        this.query = pair.Polymer0;
    //        this.reference = pair.Polymer1;
    //        this.gapScores = gapScores;
    //        makeLogLikelihood();
    //    }

    //    private void makeLogLikelihood()
    //    {
    //        loglikelihood = new double[query.Seq.Length][];
    //        logSumLikelihood = new double[query.Seq.Length];
    //        for (int qPos = 0; qPos < query.Seq.Length; qPos++)
    //        {
    //            double likelihood = 0;
    //            double sum = 0;
    //            for (int iNuc = 0; iNuc < 5; iNuc++)
    //            {
    //                sum += query.Seq[qPos][iNuc];
    //            }
    //            logSumLikelihood[qPos] = Math.Log(sum);

    //            loglikelihood[qPos] = new double[reference.Seq.Length];
    //            for (int rPos = 0; rPos < reference.Seq.Length; rPos++)
    //            {
    //                likelihood = 0;
    //                for (int iNuc = 0; iNuc < 5; iNuc++)
    //                {
    //                    likelihood += query.Seq[qPos][iNuc] * reference.Seq[rPos][iNuc];
    //                }
    //                loglikelihood[qPos][rPos] = Math.Log(likelihood);
    //            }
    //        }
    //    }

    //    override public bool IsGap(double[] n)
    //    {
    //        // this only returns true if n is a pure gap. We should consider relaxing this requirement
    //        return n[4] == 1;
    //    }

    //    override public double ComparisonScore(int posInQuery, int posInReference)
    //    {
    //        return loglikelihood[posInQuery][posInReference];
    //    }

    //    override public double OpenInsertionScore(int posInQuery, int posInReference)
    //    {
    //        return query.Seq[posInQuery][4] + gapScores.OpenInsertionScore;
    //    }

    //    override public double ContinueInsertionScore(int posInQuery, int posInReference)
    //    {
    //        return query.Seq[posInQuery][4] + gapScores.ContinueInsertionScore;
    //    }

    //    override public double OpenDeletionScore(int posInQuery, int posInReference)
    //    {
    //        return gapScores.OpenDeletionScore;
    //    }

    //    override public double ContinueDeletionScore(int posInQuery, int posInReference)
    //    {
    //        return gapScores.ContinueDeletionScore;
    //    }

    //    override public void SetQuery(Polymer<double[]> query)
    //    {
    //        this.query = query;
    //        makeLogLikelihood();
    //    }

    //    public override Polymer<double[]> GetQuery()
    //    {
    //        return query;
    //    }

    //    override public void SetReference(Polymer<double[]> reference)
    //    {
    //        this.reference = reference;
    //        makeLogLikelihood();
    //    }

    //    public override Polymer<double[]> GetReference()
    //    {
    //        return reference;
    //    }

    //    override public double DefaultQueryScore(int posInQuery)
    //    {
    //        return logSumLikelihood[posInQuery] - 1.386294;
    //    }

    //    override public double DefaultQueryScore()
    //    {
    //        return double.NaN;
    //    }

    //    override public double[] DefaultQueryScoreArray()
    //    {
    //        double[] returnValue = new double[query.Seq.Length];

    //        returnValue[0] = DefaultQueryScore(0);
    //        for (int i = 1; i < query.Seq.Length; i++)
    //        {
    //            returnValue[i] = DefaultQueryScore(i) + returnValue[i - 1];
    //        }
    //        return returnValue;
    //    }

    //    public static double[] ScoreArray(Polymer<double[]> query)
    //    {
    //        double[] loglikelihood = new double[query.Seq.Length];
    //        double[] pmfForN = new double[5] { 0.25, 0.25, 0.25, 0.25, 0 };
    //        for (int pos = 0; pos < query.Seq.Length; pos++)
    //        {
    //            double likelihood = 0;
    //            likelihood = 0;
    //            for (int iNuc = 0; iNuc < 5; iNuc++)
    //            {
    //                likelihood += pmfForN[iNuc] * query.Seq[pos][iNuc];
    //            }
    //            loglikelihood[pos] = Math.Log(likelihood);
    //        }
    //        // make cumulative
    //        for (int pos = 1; pos < query.Seq.Length; pos++)
    //        {
    //            loglikelihood[pos] += loglikelihood[pos - 1];
    //        }
    //        return loglikelihood;
    //    }

    //}

    /// <summary>
    /// Abstract class that specifies what methods alignment kernel must provide. It may compare Polymers of different nucleotide encoding types.
    /// </summary>
    /// <typeparam name="T">The type of the first nucleotide encoding.</typeparam>
    /// <typeparam name="U">The type of the second nucleotide encoding.</typeparam>
    public abstract class AlignmentKernelSymmetric<T>
    {
        /// <summary>
        /// Reports the status of the symbol as a gap.
        /// </summary>
        /// <param name="n">The symbol to be tested.</param>
        /// <returns>True if the symbol represents a gap.</returns>
        abstract public bool IsGap(T n);
        /// <summary>
        /// Computes the comparison score for the alignment of the nucleotide at a given position in the query with the 
        /// nucleotide at a given position in the reference.
        /// </summary>
        /// <param name="pos0">The position in query.</param>
        /// <param name="posInReference">The position in the reference.</param>
        /// <returns></returns>
        abstract public double ComparisonScore(int pos0, int pos1);
        /// <summary>
        /// Computes the score for continuing an insertion in the query at a given query position when that position 
        /// is aligned against a given position in the reference.
        /// </summary>
        /// <param name="posInQuery">The query position.</param>
        /// <param name="posInReference">The reference position.</param>
        /// <returns>The cost of opening an insertion at the given coordinates.</returns>
        abstract public double OpenGapScore(int pos0, int pos1);
        /// <summary>
        /// Sets the query Polymer.
        /// </summary>
        /// <param name="query">The Polymer to set as query.</param>
        abstract public double ContinueGapScore(int pos0, int pos1);
        /// <summary>
        /// Computes the score for opening a deletion in the query at a given query position when that position 
        /// is aligned against a given position in the reference.
        /// </summary>
        /// <param name="posInQuery">The query position.</param>
        /// <param name="posInReference">The reference position.</param>
        /// <returns>The cost of opening an insertion at the given coordinates.</returns>
        abstract public void SetSeq(Polymer<T> seq, int sequenceNumber);
        /// <summary>
        /// Gets the query.
        /// </summary>
        /// <returns>The query Polymer.</returns>
        abstract public Polymer<T> GetSeq(int sequenceNumber);

        abstract public double PriorScore(int sequenceNumber, int position);

        abstract public double[] GetPriorScoreArray(int sequenceNumber);
    }

    public class PMFAlignmentKernelSymmetric : AlignmentKernelSymmetric<double[]>
    {
        private Polymer<double[]>[] sequences;
        private GapScores gapScores;
        private double mutationCost = MixedNucleotideByte.LogPrior;

        public GapScores GapScores
        {
            set { gapScores = value; }
        }

        public PMFAlignmentKernelSymmetric()
        {
            this.gapScores = new LPLib.GapScores();
            sequences = new Polymer<double[]>[2];
        }

        public PMFAlignmentKernelSymmetric(Polymer<double[]> sequence0, Polymer<double[]> sequence1)
        {
            sequences = new Polymer<double[]>[2];
            sequences[0] = sequence0;
            sequences[1] = sequence1;
            this.gapScores = new GapScores();
        }

        public PMFAlignmentKernelSymmetric(PolymerPair<double[], double[]> pair)
        {
            sequences = new Polymer<double[]>[2];
            sequences[0] = pair.Polymer0;
            sequences[1] = pair.Polymer1;
            this.gapScores = new GapScores();
        }

        override public bool IsGap(double[] n)
        {
            return n[4] == 1;
        }

        override public double ComparisonScore(int position0, int position1)
        {
            return Math.Log(NucleotidePMF.IDFlux(sequences[0].Seq[position0], sequences[1].Seq[position1]));
        }

        override public double OpenGapScore(int position0, int position1)
        {
            return mutationCost + gapScores.OpenInsertionScore + PureNucleotideInt.LogPrior;
        }

        override public double ContinueGapScore(int posInQuery, int posInReference)
        {
            return gapScores.ContinueInsertionScore + PureNucleotideInt.LogPrior;
        }

        override public void SetSeq(Polymer<double[]> sequence, int sequenceNumber)
        {
            sequences[sequenceNumber] = sequence;
        }

        public override Polymer<double[]> GetSeq(int sequenceNumber)
        {
            return sequences[sequenceNumber];
        }

        override public double PriorScore(int sequenceNumber, int position)
        {
            return PureNucleotideInt.LogPrior;
        }

        override public double[] GetPriorScoreArray(int sequenceNumber)
        {
            double[] returnValue = new double[sequences[sequenceNumber].Seq.Length];

            returnValue[0] = PureNucleotideInt.LogPrior;
            for (int i = 1; i < sequences[sequenceNumber].Seq.Length; i++)
            {
                returnValue[i] = PureNucleotideInt.LogPrior + returnValue[i - 1];
            }
            return returnValue;
        }
    }

    public class MNMNAlignmentKernelSymmetric : AlignmentKernelSymmetric<byte>
    {
        private NucleotideEvolutionModel model;

        private Polymer<byte>[] sequences;
        private GapScores gapScores;
        private double distance;
        public double Distance
        {
            get { return distance; }
            set
            {
                distance = value;
                model.Time = distance;
                computeComparisonScoreArray();
                mutationCost = Math.Log(1 - Math.Exp(-distance));
            }
        }

        private double mutationCost;
        public GapScores GapScores
        {
            set { gapScores = value; }
        }

        private double[,] comparisonScore;

        public MNMNAlignmentKernelSymmetric(NucleotideEvolutionModel model)
        {
            sequences = new Polymer<byte>[2];
            this.model = model;
            this.gapScores = new GapScores();
            this.distance = model.Time;
            mutationCost = Math.Log(1 - Math.Exp(-distance));
            computeComparisonScoreArray();
        }

        public MNMNAlignmentKernelSymmetric(NucleotideEvolutionModel model, Polymer<byte> sequence0, Polymer<byte> sequence1)
        {
            this.model = model;
            sequences = new Polymer<byte>[2];
            sequences[0] = sequence0;
            sequences[1] = sequence1;
            this.gapScores = new GapScores();
            this.distance = model.Time;
            mutationCost = Math.Log(1 - Math.Exp(-distance));
            computeComparisonScoreArray();
        }

        public MNMNAlignmentKernelSymmetric(NucleotideEvolutionModel model, PolymerPair<byte, byte> pair)
        {
            this.model = model;
            sequences = new Polymer<byte>[2];
            sequences[0] = pair.Polymer0;
            sequences[1] = pair.Polymer1;
            this.gapScores = new GapScores();
            this.distance = model.Time;
        }

        private void computeComparisonScoreArray()
        {
            comparisonScore = new double[MixedNucleotideByte.NSymbols, MixedNucleotideByte.NSymbols];
            for (byte b1 = 0; b1 < MixedNucleotideByte.NSymbols; b1++)
            {
                double[] p1 = TypeConverter.ToPMF(b1);
                for (byte b2 = 0; b2 < MixedNucleotideByte.NSymbols; b2++)
                {
                    double[] p2 = TypeConverter.ToPMF(b2);
                    double x = 0;
                    for (int i1 = 0; i1 < 4; i1++)
                    {
                        for (int i2 = 0; i2 < 4; i2++)
                        {
                            x += p1[i1] * p2[i2] * model.TransitionProbability(i1, i2);
                        }
                    }
                    comparisonScore[b1, b2] += Math.Log(x);
                }
            }
        }

        override public bool IsGap(byte n)
        {
            return n == 16;
        }

        override public double ComparisonScore(int position0, int position1)
        {
            return comparisonScore[sequences[0].Seq[position0], sequences[1].Seq[position1]];
        }

        override public double ContinueGapScore(int position0, int position1)
        {
            return gapScores.ContinueInsertionScore + PureNucleotideInt.LogPrior;
        }

        override public double OpenGapScore(int position0, int position1)
        {
            return mutationCost + gapScores.OpenDeletionScore;
        }

        public override Polymer<byte> GetSeq(int sequenceNumber)
        {
            return sequences[sequenceNumber];
        }
        override public void SetSeq(Polymer<byte> sequence, int sequenceNumber)
        {
            sequences[sequenceNumber] = sequence;
        }

        override public double PriorScore(int position0, int position1)
        {
            return PureNucleotideInt.LogPrior;
        }

        override public double[] GetPriorScoreArray(int sequenceNumber)
        {
            double[] returnValue = new double[sequences[sequenceNumber].Seq.Length];

            returnValue[0] = PureNucleotideInt.LogPrior;
            for (int i = 1; i < sequences[sequenceNumber].Seq.Length; i++)
            {
                returnValue[i] = PureNucleotideInt.LogPrior + returnValue[i - 1];
            }
            return returnValue;
        }
    }



}
